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1. Introduction 

We consider Feynman integrals in D-dimensional Minkowski space with one time- and (D—l) 
Euclidean space dimensions, e =D — 4 and £ € M with |e| ^ 1, and with at most one mass. Here 
the discrete Mellin parameter n comes from local operator insertions. As worked out in detail 
in [19,25] these integrals can be transformed to integrals of the form 

^{e,n)=C{e,n,M) dxi... dx,„ 'r^ , (1.1) 

with k £'N, r\,...,rk G'N and where p{e) is given by a rational function in e, i.e., j3(£) G Q(£), 
and similarly «,• /(£,«) = ni in + cCi^i for some G {0, 1} and a,;/ G Q(£), see also [26] in the case 
no local operator insertions are present. C{e,n,M) is a factor, which depends on the dimensional 
parameter £, the integer parameter n and the mass M. Pi{x\, . . . ,x,„), Q{x\ , . . . ,x,„) are polynomials 
in the x,. Integrals of the type ( |1.1[ ) emerge in the calculation of unpolarized and polarized mas- 
sive operator matrix elements (OMEs) [2, 6, 10-14, 22, 24] and in other single scale higher loop 
calculations. In [14,22] 3-loop moments of the corresponding OMEs have been calculated. 

In addition, such integrals (1.1) can be transformed to proper hypergeometric multi-sums of 
the form' 

oo oo Li(n) L,.(n.ku--,k-\) I r(t. A r(f , ,\ 

ne,n)= I ... 11 ... I IC,(£,.,M)-iiM_^^ (1.2) 

Here the upper bounds Li{n),... ,Ly{n,ki ,^v-i) are integer linear (i.e., Unear combinations of 
the variables over the integers) in the dependent parameters or oo, and ti ^ are linear combinations 
of the n^. of the ki,... ,k^,, and of £ over Q. 

Finally, if the sums (1.2) are uniformly convergent, one of the most common tactics is as 
follows. First one expands the summand of (1.2), say 

F{n,ni,...,nr,vi,...,Vk) = F,{n,ni, . . . ,Vk)e' + Ft{n,ni,. . . ,Vk)e'^^ + . . . (1.3) 

with f G Z by using formulas such as 

r(n + ! + £) = ^^"^^1^ t^f and Bin, 1 + £) = 1 exp 

B{n,l + e) n y^^j 

(1.4) 

with £ = r£ for some r G Q. Here B{x,y) = r{x)r(y)/r{x +y) denotes the Beta-function and 
•^1, . . . ,1 («) is a special instance of the harmonic sums [15, 40] defined by 

" sign(mi)'i sign(m^)'* 

(1 = 1 ij 1^=1 Ij^ 

with mi,...,mk being nonzero integers. Then one applies the summation signs to each of the 
coefficients in (1.3). I.e., the ith coefficient of the £-expansion of ( [L^ ) yields 

oo oo L[{n) Ly(n,k\,...,ky-\) I 

L ••• L L ••• L ^/^•(«,«i,...,«„vi,...,v^). 

n\ = \ n,-=lk\ = i ky=\ k=l 




'For convenience, we assume that the summand is written terms of the Gamma function r(x). Later, also Pochham- 
mer symbols or binomial coefficients are used which can (if necessary) be rewritten in terms of Gamma-functions. 
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Then the essential problem is the simplification of these sums to special functions, like, e.g., har- 
monic sums, 5-sums [30] 



n 1 1 Ik- 1 y'k 

^mi,...,mi- {^It ■ ■ ■ j^k^f^) — 2^ "Imf • • • ^ "jm^ i (1-6) 
ii=l h ij,=l h 

cyclotomic harmonic sums [4], or more generally to indefinite nested sums and products [38]. For 
various special cases, this simplification can be carried out with efficient methods available, e.g., in 
Form; see in [30,40]. 

More general sums can be handled with the Mathematica package EvaluateMultiSum [3,23] 
based on the summation package Sigma [36]. With the underlying difference field algorithms [7, 
29,34,35,37,38] generalizing the hypergeometric summation paradigms [32] to multi-summation 
we are currently simplifying sums up to nesting depth 7. The compact representation of the output, 
i.e., the elimination of all algebraic relations among the arising indefinite nested sums and products 
can be guaranteed by difference field theory [39]. For harmonic sums, cyclotomic sums, S-sums 
and cyclotomic 5-sums and their infinite versions (quasi-)shuffle algebras are utilized; see e.g., [4, 
5, 17,20,21] and references therein. In this regard, the Mathematica package HarmonicSums is 
heavily used [1]. This general machinery has been applied to non-trivial massive 3-loop diagrams 
arising, e.g., in [2,6, 12,24]. 

Another possibility is the method of hyperlogarithms [27] which can be used to evaluate integrals 
of the form ( |l.l[ ) for specific values ?i G N if one can set e = 0. An adaption of this method for 
symbolic n has been described and applied to massive 3-loop ladder graphs in [2]. 

In this article we follow another promising approach for the all n expansion 
1. Calculate a recurrence in n for the multi-integral ( |1.1[ ) or multi-sum (1.2) 



2. Given this recurrence and initial values (using, e.g., the EvaluateMultiSum package), 
calculate the £-expansion by a recurrence solver for Laurent expansions. 

In [25] we followed this approach by calculating recurrences of multi-sums using techniques of [42] 
and efficient algorithms developed in [41]. Subsequently, we present two new techniques to com- 
pute such recurrences. In the first approach we apply an enhanced and optimized version [1] 
of the multivariate Almkvist-Zeilberger algorithm [9] to calculate recurrences for Feynman in- 
tegrals in the form ( |l.l[ ). In the second approach we use a common framework [33] within the 
summation package Sigma that combines difference field [29, 35] and holonomic summation tech- 



niques [28,43] to compute recurrences for Feynman integrals in the form ( |1.2[ ). Two new packages, 
Multi Integrate by J. Ablinger and RhoSumby M. Round facilitate these tasks completely 
automatically for the integral or sum representation, respectively. 

The outline of the article is as follows. In Section ^ we present the general idea of how the 
Laurent series representation can be calculated from the given recuiTcnce representation. With this 
knowledge, we illustrate our integration and summation methods in Sections ^ and 0, respectively. 



2. Finding Laurent series solutions of linear recurrences 

One of the key ingredients of the summation and integration tools under consideration is a 
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recurrence solver for £-expansions. To illustrate the ideas of the solver, we consider the single sum 

+k,n) _ - r(n)r(f +^) j_ ^ ^ ^ ^ 



■.Fo{n)+Fi{n)e+F2{n)e^ + ... (2.1) 



which is related, e.g., to sums arising in [12]. In order to derive the coefficients Fi{n), we compute 
in a first step a recuiTcnce relation using the summation package Sigma. Internally, it activates a 
difference field version of Zeilberger's creative telescoping paradigm [32]. In our example it turns 
out that o5^(e,?i) satisfies for all integer n > 1, as an analytic function in e throughout an annular 
region centered by 0, the recun^ence 



ao{e,n)y{e,n)+ai{e,n)y{e,n+l)+a2{e,n)y{e,n + 2) 

_ -4{2n + 3) (£ + 4«2 + 12« + 8) r (f + l) r{n+l) 
~ (n+l)(« + 2)2(e + 2n + 2)r(f + « + !) 

with ao(£,«) = —4n{n+ 1), ai{e,n) = 4{n+ l)(e — 2« — 3), and a2{e,n) = (e - 2?i — 4)^. 
we compute the e-expansion /io(n) +/ii(«)e + h2{n)e^ + ... of the right hand side of (| 
formulas such as ( L4); here we get /io(«) 



(2.2) 

Next, 
by 



(2n+3)^Sl(n) , (2n+3)Si(nf (2n+3)52(«) 



(«+l)2(n+2)2+ (n+l)(«+2) {n+l){n+2) 



+ 



(«+l)(«+2)'"U"; 
(2n+3)2 



f 4{2n+3)S,{n) 



(»+l)(«+2) +(n+l)2(n+2) 



2(2»+3)^ 



-),h2{f 



+ („^i)3(„^2)^ ) • '^^ ^ consequence, for the Laurent series 



expansion 5^{E,n) = Fo{n) +Fi{n)£ + .. . the following relation holds: 

ao{£,n) FQ{n)+Fi{n)e + ...]+ai{£,n) Fo{n + \) + Fi{n+ \)£ + . . .] 
+a2{£,n) Fo{n + 2)+ Fiin + 2)£ + ...]= ho{n) + hi{n)e + h2{n)e^ + ... 
Next, we expand the first two initial values^ of (£,«),« = 1,2: 



(2.3) 



2-C2 + (f-C2)£ + (- 



-^ + | + 1)£2 + . 



^(s,2) = |-| + (^-^)s + (H^-:i|-^)s^ + . 



(2.4) 



' 4 32^^ + 1 32 16 

by using the package EvaluateMultiSum [23] in Mathematica; alternatively, one could use the 
package Summer [40] in Form. Then given the recurrence (2.3) and the first initial values (2.4) 
(to be more precise the polynomials ai{e,n), the first coefficients hi{n), and the first values of their 
expansions), we are ready to calculate the first three coefficients of the all n series expansion (p|). 
Namely, by setting £ = in (23), it follows that the constant term Fo(?i) satisfies the recurrence 

ao(0,«)^b(«)+ai(0,?i)fo(«+l) + a2(0,«)fo(« + 2) = /zo(«). (2.5) 

Note that together with fb(l) =2 — 1^2 and Fq{2) = ^ — ^ the sequence Fo(n) is completely deter- 
mined. At this point we exploit algorithms from [7, 8, 31, 35] which can constructively decide if a 
solution with certain initial values is expressible in terms of indefinite nested products and sums. 
More precisely, with Sigma one obtains 



(2.6) 



^Here stands for the Riemann Zeta function i^(r) = 
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Now, plugging in the partial solution 



into ( [2.3| ) and moving F(){n) to the right hand side yields 

ao{e,n) [Fi {n)£+F2{n)E^ + ...]+ (£,«) [fj {n + 1)£ + F2{n +1)e^ + ...] 

+ a2ie,n)[Fi{n + 2)e+F2{n + 2)e^ + ...] = hQ{n)e + h[{n)e^ + . . . (2.7) 

with h' (n) - 4(2»+3)Si(«) , 2{4n+5) , 
With ho[n) - („+i)(„+2) + {n+\)Hn+2y 

_ -(2»^+9»+8)(2»+3) _ (2„+3)^Si(«) _ {2n+3)Si{nf (-2«-3)52(«) , (-l)"(25-2(«)+C2) 
"1'^"^ («+l)^(«+2)3 («+l)2(«+2)2 (n+l)(«+2) (n+l)(n+2) «+2 



Observe that the coefficients and the inhomogeneous side of the recurrence (2.7) have all the factor 
£. Hence dividing the recurrence through e, we obtain again a recurrence of the form (2.3) where 
the explicitly given hi are changed to /j- and Fi is the new constant term. In other words, we can 
repeat this construction process for Fi{n). Namely, by coefficient comparison Fi{n) is uniquely 
determined by 

ao(0, n)Fi (n) + ai (0, n)Fi («+!)+ ^2(0, n)Fi {n + 2)= h'^in) 



and the initial values given in ([2.4|). Solving the recurrence with these initial values leads for all 
n > 1 to the sum representation 

FJn) = (-lyU^-^^ - !^zMM) + (-l)"C2^iW + 2{-irSdn)S^2{n) _ ^2 8) 

2n n ^ n n 

Similai^ly, one can loop further and calculates, e.g., the coefficients FQ{n),. . . ,F(,{n) (in terms of (^2> 

(^3, (^5, (^7 and harmonic sums up to weight 8) in about 30 minutes. 

More generally, let o5^(£,n) be a function 

• which is for each integer n with n> X analytic in £ throughout an annular region centered 
by 0; let ^{e,n) = TT=t^iin)£' be its Laurent series expansion for some ? € Z; 

• which satisfies (as £-expansion) the recurrence 

ao{e,n)y{e,n) H \-ad{e,n)y{e,n + d) = /i,(?i)£' +/if+i(«)£'+^ ... 

for polynomials ai{e,n) eK[e,t] with a^/ (0, n) /Ofor alln > A and for functions where 
the first coefficients ht{n),hf^i (n) ... ,hu{n) with « > A are expressible in terms of indefinite 
nested sums and products. 

Then there is the following algorithm [25, Cor. 1] implemented in Sigma: 

Input: the polynomials a,(£,«) and the product-sum expressions hj{n) {t < i < u) as above; the 

values Cij {t <i <u, X < j < X + d) such that y{e,j) = ctjs' + Ct+ijs'^^ H h c„,y£" + . . . 

Output: the maximal r G {—oo^t,t + l,...,u} such that the coefficients Ff(?i),F,+i (n) ... of 
the £-expansion of y{e,n) can be expressed in terms of indefinite nested sums and products. If 
r / —00, these representations of the coefficients are computed explicitly. 

For rigorous proofs, further details concerning efficiency, generalizations, and the function call 
within the package Sigma we refer to [25]. 
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3. Calculating e-expansions for multi-integrals 



We aim at computing a recurrence relation for multi-integrals of the form ( |1 . 1| ) and finding a 
Laurent-series solution that agrees with the input integral. Here we exploit an enhanced version [1] 
of the multivariate Almkvist-Zeilberger algorithm [9] which contains as input class these integrals. 
More generally, it can handle integrands being hyperexponential in the integration variables x, 
(i.e., the logarithmic derivative of the integrand w.r.t. Xj is a rational function in the x, and n) and 
hypergeometric in the discrete parameter n (i.e., the shift quotient w.r.t. n of the integrand is a 
rational function in the x, and n). 

In order to illustrate the basic ideas, consider the double integral 

J^{e,n)= f JxiJx2=fo(n)+fi(n)£+f2(«)£^ + -.- (3.1) 

■w Jo (1+Xij'^ 

^ V ' 

F{n,Xi^2)'-= 

First, one applies the multivariate Almkvist-Zeilberger algorithm. To be more precise, given d G'N 
one looks for polynomials and rational functions Ri{n,x\ ,X2) such that 

eQ{n)F {n,x\,X2) + e2{n)F {n+l,x\,X2) H \- ed{n)F{n + d,xi,X2) 

= D.vi (/? 1 F («, xi , X2 ) ) + D;t2 (^2^ («, -^1 , -^2 ) ) ; 

here D^. stands for the differentiation w.r.t. x,. Internally, a clever ansatz is performed with undeter- 
mined coefficients which amounts to solving a linear system of equations. To hunt for a solution, 
one starts with d = and increases the recurrence order d step by step until a solution is found. In 
our particular case, for J = 1, one gets 



{n + l)F{n,xi ,X2) + {n + 2)F{n + l,xi ,X2) = D,, + D,, ^^^^J^^^ll^ . (3.2) 



Applying now the two integral signs on both sides of ( |3.2[ ) leads to the following recurrence 



-(« + 1)J^(£, «) + (« + 2) J^(£, «+ 1) = [\xi + l)"+^-''dxi- [\dxi. (3.3) 

Jo Jo 



hie,n) 



Note that the right hand side of (3.2) consists again of an integral. However, the nested depth is 
decreased from two to one. By recursion, we treat now this simpler integral again by the method 
under consideration. In this case, we find 

Wo ^_ 2"+^-l , (-2"+^(log(2)(«+2)-l)-l) ^ , (2"+'(log(2)^(«+2)^-21og(2)(«+2)+2)-l) ^2 , 
h[£,n) - + £ + ^j;^ £ +... 

Finally, plugging this result into ( ^j| ) gives a recurrence that fits into the input class of the re- 
currence solver presented in Section ^. Together with the expanded initial values (which can be 
calculated easily) 

^(£,0) = f + (I - 21og(2)) £ + (11 + logHl) - ,2 ^ _ , _ 
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we ai^e in the position to calculate the first thi^ee coefficients Fi{n) of the expansion (3.1). 

In order to calculate this integral within Mathematica, several packages have to be load in: 
the Sigma package to use the recurrence solver from Section^ and the EvaluateMultiSum 
package to deal with expansions. Finally, we load in the package Multi Integrate which 
contains an efficient implementation of the multivariate Almkquist-Zeilberger algorithm with the 
help of homomorphic image testing; see [1]. In addition, it combines all the steps described above 
(together with variations of the presented method) to perform the £-expansion. 

in[ii:= << Sigma. m 

Sigma - A summation package by Carsten Sclineider © RISC 

in[2]:= << EvaluateMultiSums.m 

EvaluateMultiSums by Carsten Schneider - © RISC 

In [3]:= << Multilntegrate.m 

Multilntegrate by Jakob Ablinger - © RISC 

Now we are ready to calculate the expansion of the integral above. The coefficients are returned in 
list form, i.e., {FQ{n),Fi{n),F2{n)}: 



in[4i:= sol=mAZExpandedIntegrate[ ^^,'^^' , n, {e, 0, 2}, {{xi, 0, 1}, fe, 0, 1}}] 



0Ut[4]: 

' (n+ 1 
^ ' ln2' 



Note that the involved sums can be rewritten in terms of S-sums, see ([E^, using the command 

Transf ormToSSums from the package Harmonic Sums: 

in[5]:= << HarmonicSums.m 

HarmonicSums by Jakob Ablinger - © RISC 

in[6]:= TransformToSSums[sol] 

.Si(2,n) Si(n) , 2''+' 1 , /-Si(2,n) 2°+^ ^ S2(2,n) SjH , (2°+l - 
Out[6]= \ + — ^ - — , ln2 \ I V ' V , \ 



n+1 n+1 (n+l)2 (n+l)^' V n+1 (n+l)^/ n+1 n+1 (n+l)^ 
^2(^l(^,^_\,, S2(2,n) 2"+l X S3(2,n) S^jn) 2°+! 1 ^ 

M2(n+1) (n+l)2y I n+1 (n+l)^/ n+1 n+1 (n+l)* (n+l)*^ 

The proposed method extends to integrals with higher nesting depth. We conclude this approach 
by the calculation the coefficients Fo{n),Fi{n),F2{n) of the £-expansion of the following triple 
integral: 

/' t d^^dx2dx3=Foin) + Fy{n)e + F2{n)e^ + ... 

Jo Jo Jo {l+X3)^ 

Integrals of this type emerge as partial integrals, e.g., for 3-loop ladder topologies. 

in[7]:= mAZExpandedlntegrate [ /".'^"'^^ ^ , n, {e, 0, 2}, 

{{xi , 0, 1}, {x2 ,0,1}, {x3 ,0,1}}]/ /TransformToSSums//ReduceToBasis 

n,r7, r -3n-4 , 2-'+H3n+4) Sjn) Si(2.n) 2-+^ (3n^+6n+2) 4n^+9n+4 , 2Si(n) 

t (n+l)^(n+2)^ + (n+l)^(n+2)^ (n+l)(ii+2) (n+l)(n+2) ' (n+l)»(n+2)^ + (n+l)^(n+2)^ (n+l)(n+2)^ 

S2(n) _ (-l)°S-i(n)+2Si(2,n) , _ 2°+H3n+4) _ Si(2,n) (-1)° , 1 N , 82(2,11) 

(n+l)(n+2) (n+l)(n+2)= "t""!^!, (n+l)2(n+2)2 (n+l)(n+2) (n+l)(n+2)= ^ (n+l)(n+2)2 J (n+l)(n+2) ' 

/ 2°(3n+4) Si(2 n) (-1)° 1 g2 , / -3n-5 2°+' (311^+611+2) g^(n) 

V (n+l)2(ii+2)2 2(n+l)(n+2) 2(n+l)(n+2)2 2(n+l)(n+2)= ^ "'^ \ (n+l)2(n+2)3 (n+l)='(n+2)3 (n+l)(n+2)2 ~^ 

-W^f ^-'W I 1 1 I 2Si(2.n) S,(2,n) n. ^ Si(n)^ (3n+4) (3n^+9n+7) 
+ l(n+l)(n+2)i + (n+l)(n+2)3j + (n+l)(n+2)2 " (n+l)(n+2) " 2(n+l)(n+2)i (n+l)4(n+2)4 
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2°+i(7n3+30n2+44n+22) 2(2n+3)Si(n) S;(n)-S2(2,n) 83(11) 2(2n+3)Si(2.n) 83(2,11) 

+ (ll+l)«(ll+2)« (n+l)2(ll+2)3 + 2(n+l)(ll+2)2 (ll+l)(ll+2) + (ll+l) = (ll+2)3 + (ii+l)(ll+2) 

81.1(1.2,11) , ^ ,^n^ 8-i(n) , 2S_2(ii) , S_i,i(ii) S2(-2.ii) 8i.i(-1.2,n) 

+ (ii+l)(n+2)2 \ (ii+l)(n+2)3 + (n+l)(n+2)2 (n+l)(ii+2)2 (n+l)(n+2)2 (n+l)(ii+2)2 J / 

Here the command ReduceToBasis ensures that no algebraic relations exist among the arising 
S'-sums; see also [1,4,5, 17,20]. 



4. Calculating £-expansions for multi-sums 

In our second strategy we rely on a common summation framework [33] of the difference field 
and holonomic approach [9,28]. The input class of the proposed method covers multi-sums of the 



form (1.2). More generally, the summand itself can be an expression in terms of indefinite nested 
sums and products. To illustrate this approach we aim at calculating the first coefficients of the 
double sum 

"^"iji^^ {-\y(-''+''--)r{j+k+i]r(n-k)(i-^) (2-f) , , 

ito pi r(-<:+„-l)(3-e),+,(f+3]^.^^. 

' V ' 

F(n,k) 

First one considers the inner sum f («,^) which itself is an analytic function in e for each integer n 
and k with n>3 and < ^ < ?i — 3. First, we hunt for a recurrence relation in k. As described in 
Section § one can use the package Sigma and calculates 

{k+\){e-2k-2){k-n + 2)F{n,k) 

- {k-n+\){e^ -ek + en-e-Ak^ + 2kn- Uk- U)F{n,k+ \) 
-2{e + k + 2){k-n + \){k-n + 2)F{n,k + 2) = eo{n,k) + ei{n,k)e + e2in,k)e^ + . . . (4.1) 

with 

l6{k-n + \){k-n + 2){k^n + 4k^ + 2kn + 22k-4n + 2S,) 



eo{n,k) = - 
e\ {n,k) = 



(k + 2)^k + 3)^k + 4f 

{k-n+i){k-n + 2) {2k^n+ lOk^ + IS/fn + n4k'*+43k^n + 69&k^ - 32k^n + n60k^ - 220kn + 2l40k - 184« + 1000) 



(^ + 2)3(i + 3)3(yt + 4)3 

, , 1, , , ^,\6(k^n + 4k^+2kn + 22k-4n + 2%)S2{k) 

e,in,k) = --i^k+mik-n+m-n + 2)i ,^^,^,^.^,^,y^,^,yl 

-2{5k^°n + 4k^° + 142/t'« + 230/t' + 1602*:^i + 3812/1:'* +9552/t''n + 31332/t'' +32861/t*n 

+ 150820*:'* + 64970/t^« + 452886A:^ + 64124k* n + 851364k* + 1 1424*:3n + 985264A:' - 29328/t^« 
+ 606320/t2 - 1171 2/t«+ 121344*: + 5952«- 30144) ^ {{k + \f {k + 2)\k + 3)* {k + 4f k\)) . 

In addition, one computes a mixed recun^ence, i.e., besides shifts in k we allow in addition one shift 
in n, but keep k unchanged. Using again Sigma, one obtains 

{k -n)[-e^ + ek-en + 21^ - 2kn + 2k + Ir? + 2n + 2)F{n,k) 

+ {e - n - l){e + 2n + 2){k - n + l)F{n + l,k) 
-2{e + k+l){k-n){k-n+l)F{n,k+l)=fo{n,k)+fi{n,k)e+ f2{n,k)e^ + ... (4.2) 
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with 



Mn,k) . 



I6{k-n){k-n+l){k^ -2kn + 2k-4n-i) 



Mn,k) = — 



{k+l)(k + 2)2{k + 3)^ 
-n){k-n+ 1) (2/t5 - Ak'^n + 120 - 26Pn + \7k^- 56k^n -llk^- A2kn - 47k - An - 19) 
(k + T)2(k + WW+W 
16(k^ -2kn + 2k-4n-\)S2(k) 



Mn,k) = --ik^rMk-n)ik-n + l){- ,^,^,^,^2^,^,>,,^ 

+ 2(5k^ - U)k'n+\02k' - 188<:*« + 745A:^ - 1300/t\ + 2644*:' - + 4855*:" - 8162A-^i 

((*:+l)''(<: + 2)^(*: + 3)'*A:!)). 



We emphasize that these two recun^ences together with the initial values F(4,0) = — |f + §f — 
. . . and F{4,1) = —j^ + ^ — + • • • enables one to calculate the first three coefficients of 
the e-expansion for each F{n,k) with n>4 and < ^ < « — 3. 

Given the two recurrences above (with this particular shape of shifts), one can apply the algo- 
rithm from [33] to calculate a recurrence for y{s,n). Using again Sigma, one gets the relation 



£(£ + 2)2(2n + 1)(£ - n - 3)(£ - n - 2)(£ + 2n + 2){e-n- 1) V(£,« + 1) 

+ 4£(£ + 2)(?I+l)^(?l + 3)(2« + 3)(£-?l-2)(£-«- l)^(£,?l) 
= hQ{n) + /ji {n)e + h2{n)e^ + {n)F (?i,0) 

+ hi,{n)F («, 1 ) + /J5 {n)F («, « - 3) + /i6 {n)F {n,n-2) + ... (4.3) 

with 



ho{n) 
hi{n) 
h2{n) 



64 ( I5n> + 5;!** - 188n'' - 194n'' + 655n' + 1 129;!-' - 266;!^ - 1660n2 - 1224« - 288) 
■ 9(«- ' 

32(l5n'-'+87n'--45n"-329n"'-215n''+249n*+2331n'+1545n''-6396n5-5440n''+6686/i3+7776n2+648H-864) 
" ~ 27(n-l)-'n-' ' 

1 6 (4« ' 2 + 1 6,1 1 1 - 23« _ 1 ssn' - 2 1 9«* + 208n'' + 503n* + 1 2«5 - 445«-* - 9(m^ + MOn- + 272n + 64) 52 (h ) 



(«-l)2«2 

- (^4(708,j'^ + 916n'-'^ -9426n'''- 15711«'^V30903n'- + 61829H" -49101,!'° - 115105,!' +44049n* 
+ 121347,1^- 789,1^ - 57380,i' -41832«''- 14040«^ + 16848«^+ 6048,1 -3456)) / (27(,i- 1 )''«''), 

, , , , 4E2(2n''+3«5-14,!-'-9/r'-49n--I15n-62)(H+l)- 8E(2n'' + 10/r'-15ii--30n+8)(n + l)-' 32 (3/r + lOn+s) (,!+ 1 )\ 

h3{n)=( ^ '- + ^ —I '- ^ ^ ). 

/i4 („)=- 32(3,1 + 4)(,i+1)^+8(2m^- 6,1' -25,1- 12) (,7 + l)^e- 4(2,7^ -3,!'* + 16/1^ +32,1' -18,!- 31) (,7+l)^e-, 
ftj (,i) = 1 6(n - 2) 2 (3,1^ + ,i2 - 7m - 4) (,! + 1 ) ' - 4 (2,,' - 1 5,!* + 46,7^ - 1 8,1'' - 1 25,!^'' + 86,!^ + 96,7 - 1 6) (,! + 1 )^ e 

+ 2 ( len** - 58,1^ - 7,1^ + 22 1«^ - 59,!'* - 325,i^ - 26n- + 146n + 60) , 
hein) =32(3«'' - 8n' - lOn^ +20,;+ 16) ('7+1)^ - 8(2,,'^ - 13,7^ +27n'' +56,!^ + 14,i^ + 66n + 60) {n+\fe 

+ 8 (6,1'' - 23,1* - 43«^^ + 46,1'* + 68«3 - 23,7^ - 43,i - 1 2) ; 



for the different command calls within Sigma we refer to [33, 36]. As indicated above, the algo- 
rithm itself only uses the two recun^ence relations ( 4. 1 ) and ( p^ ) and thus the occurring expres- 
sions F{n,0),F{n, l),F{n,n — 3),F{n,n — 2) remain unevaluated. Next, one applies the proposed 
method recursive on these sums. Since these objects are simpler than the input sum J>^{e,n), the 
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termination of our method is guaranteed. E.g., the calculation of the e-expansion 

._'i^ (-im;+i)rw(";^)(2-f)^ 
.il r(„-i)(3-e),(f+3)^. 

n-l( S,{nf 3Si(n)^ ( 15 352(») N c / n , 352(«) 253(») 1252.l(») ^^2 



(n— 1)« 



boils down to the method described in Section Similarly one proceeds for F{n,l) and F{n,n- 
3),F{n,n — 2). This finally leads to the following simplified right and side of (|4.3|): 



Ox£"+ (^32n(n + l)2(3«+4)(«2 + 3n+4)-128(« + l)^(3n+4)5i(«))e 

+ ( - 32(2«2 + 7n - 2)(« + («) + 16(4«^ + 12«'* - 27«^ - 129«' - 130« - 28)(n + 1)^52(«) 

- 8n (12/?*' + 58n^ - 27n'* - 566n^ - 1 125n2 - 804n - 156) (« + 1 ) - 32(3/1 + 4) (n + 1 £2 - 



Note that the left hand side of ( p3| ) and its right hand side (after the simplification) can be divided 
by £, i.e., the coefficient of J>^{e,n +1) evaluated at e = does not vanish. Hence together with the 
expanded initial values y{e,3) = — | + jjS + . . . and ^(£,4) = + ^£ + • • • one can activate 
our recurrence solver for Laurent series to calculate the first two coefficients of the £-expansion 

cfifr. \ 8 c ^ \ /l^l , i\c ^ \ , 4)i(3«+7) , /-2(2«3+3n2+3n+6) , , 

(6n3+23n2+27n+12) 2 c ^ ^2 o , 1 \ c / \ «(«+l)(5«+14) \ ^ , 

+ („+i)(„+2) '^2(^) - („+i)(„+2) '^l (") - 2(2" + 1)^3 («) - („^2)^ j ^ + • • ■ 

Summarizing, in the presented method one constructs step by step suitable inhomogeneous recur- 
rences from the innermost sum to the outermost sum^. As one can see already for double sums, this 
construction is quite involved and is fairly complicated for more nested sums (e.g., taking care of 
poles, estimating how far one should expand^, or exploiting a refined difference field theory [37]). 
The new package RhoSum deals with all these aspects using as backbone the packages Sigma, 
Harmonic Sums, and EvaluateMultiSums. After loading 

in[8]:= << RhoSum.m 

RhoSum - Package for Refined Holonomic Summation © RISC 

we can perform the calculation from above with the function call 

,n[9i:. FindSum[ ^ ' J '^"^ , {{j, 1, n - k - 2}, {k, 0, n - 3}}, {n}, {3}, 

r(-K+„-i)(3-.w(,+3)^^^ Expanding {ep,0,l} 

{-CTT2)Sl(n)-4(2n+l)S2(n) + 5H^, 



(n+l)(n+2)' 

-2(2n^+3n^+3n+6j (6n3+23n^+27n+12j 2 / ^,2 oCo^^l^q (..\ n(n+l)(5n+14) ^ 
(n+l)^(n+2)^ Si W + (n+l)(n+2) - J^^^Tji^)^^'''') -2(2n+ 1)83(11) ) 



^^In [28] this idea has been considered for homogeneous recurrences with polynomial coefficients. In our ap- 
proach [33] we observed that setting up the recurrence system in the special form given above (instead of allowing 
a general holonomic system) one can derive an efficient algorithm without using Grobner basis. In this way, the holo- 
nomic approach could be extended in [33] to handle also inhomogeneous recurrences formulated in difference fields. In 
order to take into account the e-expansion of the inhomogeneous sides, new ideas have been added into Sigma. 

"^E.g., in the illustrated example from above one has to start to calculate three coefficients of the e-expansion and 
ends up only with the first two coefficients. 
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A more involved problem is, e.g., 



"_2 j n-4+j,-j (_i)-;.-;2+«-v¥(^g(-^+^.--4)r(-|-;+;,-;,+„^ ^ 
jfe) jhl (e+l)(e+2)r(2-f)r(e+4)r(y-7,+72+2)r(-2e-y+7,+n) ^ 

r(f+3)V(-^)r(-f)r(e)r(,Wi+i)r(-g+7i+i)r02+i) i 3,^ ^-2,^ p-i,^ , 
r ( f +« j r ( - f -j+ji -72 +« j 

Sums of this type occur, e.g., in case of 3-loop topologies with one massive and one massless 
fermion line. If we insert this sum into Mathematica in the variable /, then we get the following 
expansion. The constant term is too large to present it here. 

in[ioi:= FindSum[f,{{j2,l,n-4+ji-j},{ji,0,j},{j,0,n-2}},{n},{2},M,ExpandIn^{ep,-3,-l} 

n„trin, f 4(n='-5n^+6n-4) 2(-l)°(n^-2n3+n^-12n+8) 8(-l)°(n-3)Si(n) 
""tnoi- t 9(n-2)2(n-l)n3 + 9(n-2)(n-l)n3 9(n-2)n2 

16(-irs_,(n) , (-l)°(10-n)Si(n)^ , . (-l)" (l7n^-179n^+240n-96) 4 . 

^ 9(n-2)n2 ^ 9(n-2)n= \ 27(n-2)(n-l)n3 ^ 9(n-2)2(n-l)nV ^1 W 

-l)''(-2nHl9n^-5n^ + llln=-190n+84) ^ 2 (8n^-63n* + 163n'-218n^ + 140n-36) ^ (-l)°(14-n)S2(n) 



27(n-2)(n-l)n« ^ 27(n-2) = (n-l)=n* 

(-l)°(22-n)Si(n)^ (-l)-(n3-299n^+412n-168) 1 N., ,.2 

54(n-2)n2 ^\ 108(n-2)(n-l)n3 ^ 9(n-2) = (n-l)nV ^1 '-'v 

/ -6nH5n-5 , (-1)° (-313n^+1708n3-2262n^ + 1740n-648) (-l)»(7n+46)S,(n) N „ , ^ 
l27(n-2) = (n-l)2n2 ^ 162(n-2)(n-l)n4 ^ 18(n-2)n= ^^IW 

I -289n'^+2430ll^-6803n''+10222n^-8536n^+4176n-1008 
^ 324(n-2)2(n-l)=n= 

(-l)»(475n^-1768n'^+2423n=-6170n''+17912n^-21856n^ + 15744ll-5184) 2(-l)°(lln-18)S2 i(n 



^ 648(n-2)=(n-l)n5 9(n-2)n2 

,„ ^„^^ 16(-l)°Si(n) 4(-l)°(9n^+22n^-43n+18) ^ 2(-l)°(3n-16)S_3(n) 
+ 0-2Wl^ 9(n-2)n 27(n-2)(n-l)na 9(n-2)n^ 

^ -6n^+30n^-35n+24 , (-1)° (n^ -641n3+1646n^-976n+192) ^ ^ , (-l)°(47n-8)S3(n) 

+ 9(n-2)2(n-l)n3 ^ 108(n-2)2(n-l)n3 J^2Vnj+ 27(n-2)n2 

4(-l)°S_;.i(n) ^ -n^+5n=-6n+4 , (-1)° (n'-2n^+n^-12n+8) ^ (-l)°(3-n)Si(n) ^ ^ \ 
3(n-2)n + U(n-2)2(n-l)n3 ^ 12(n-2)(n-l)n^ ^ 3(n-2)n^ J ^2 J 

The above expressions can be analytically continued to complex values of the Mellin variable n 
using relations given in [16, 18-20]. 



5. Conclusion 

Massive Feynman integrals with operation insertion can be expressed in terms of multi-in- 
tegrals and multi-sums over hypergeometric and hyperexponential functions. We presented new 
methods to calculate the first coefficients of the £-expansion of such multi-sums and multi-integrals. 
Here the multivariate Almkvist-Zeilberger algorithm and the common framework of the holo- 
nomic and difference field algorithms have been enhanced to calculate recuiTcnces. Then a re- 
currence solver for Laurent series expansion is used to extract the all n coefficients of the e- 
expansion. Besides of the usage of the Mathematica packages Sigma, Harmonic Sums and 
EvaluateMultiSums, two new packages Multilntegrate and Rho have been developed 
that can carry out these calculations in a completely automatic fashion. 

Acknowledgment. This work has been supported in part by DFG Sonderforschungsbereich Tran- 
sregio 9, Computergestiitzte Theoretische Teilchenphysik, Austrian Science Fund (FWF) grant 
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